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r—i Abstract 

We describe a method of asymptotic approximations to solutions of mixed boundary value problems for 
the Laplacian in a three-dimensional domain with many perforations of arbitrary shape, with the Neumann 
boundary conditions being prescribed on the surfaces of small voids. The only assumption made on the 
geometry is that the diameter of a void is assumed to be smaller compared to the distance to the nearest 
neighbour. The asymptotic approximation, obtained here, involves a linear combination of dipole fields 
constructed for individual voids, with the coefficients, which are determined by solving a linear algebraic 
system. We prove the solvability of this system and derive an estimate for its solution. The energy estimate 
is obtained for the remainder term of the asymptotic approximation. 
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1 Introduction 
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In the present paper we discuss a method for asymptotic approximations to solutions of the mixed problems for 
the Poisson equation for domains containing a number, possibly large, of small perforations of arbitrary shape. 
The Dirichlet condition is set on the exterior boundary of the perforated body, and the Neumann conditions are 
specified on the boundaries of small holes. Neither periodicity nor even local "almost" periodicity constraints are 
imposed on the position of holes, which makes the homogenization methodologies not applicable (cf. Chapter 
4 in [5], and Chapter 5 in 9 ). Two geometrical parameters, e and d, are introduced to characterize the 
maximum diameter of perforations within the array and the minimum distance between the voids, respectively. 
Subject to the mesoscale constraint, e < const d, the asymptotic approximation to the solution of the mixed 
boundary value problem is constructed. The approximate solution involves a linear combination of dipole fields 
constructed for individual voids, with the coefficients determined from a linear algebraic system. The formal 

>— : asymptotic representation is accompanied by the energy estimate of the remainder term. The general idea of 

mesoscale approximations originated a couple of years ago in [5J , where the Dirichlet problem was considered 

y—i for a domain with multiple inclusions. 

The asymptotic methods, presented here and in [6], can be applied to modelling of dilute composites in 

• ^h problems of mechanics, electromagnetism, heat conduction and phase transition. In such models, the boundary 

conditions have to be satisfied across a large array of small voids, which is the situation fully served by our 
approach. Being used in the case of a dilute array of small spherical particles, the method also includes 
the physical models of many point interactions treated previously in [H [2j E] and elsewhere. Asymptotic 
approximations applied to solutions of boundary value problems of mixed type in domains containing many 
small spherical inclusions were considered in [I] . The point interaction approximations to solutions of diffusion 
problems in domains with many small spherical holes were analysed in [5] . Modelling of multi-particle interaction 
in problems of phase transition was considered in [3] where the evolution of a large number of small spherical 
particles embedded into an ambient medium takes place during the last stage of phase transformation; such 
a phenomenon where particles in a melt are subjected to growth is referred to as Ostwald ripening. For the 
numerical treatment of models involving large number N of spherical particles, the fast multipole method, of 
order O(N), was proposed in [3J, and it appears to be efficient for the rapid evaluation of potential and force 
fields for systems of a large number of particles interacting with each other via the Coulomb law. 
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We give an outline of the paper. The notation fl^ will be used for a domain containing small voids F"', 
j = 1, ... , N, while the unperturbed domain, without any holes, is denoted by f2. The number N is assumed 
to be large. 

If Q is a bounded domain in M 3 , we introduce L 1,2 (£l) as the space of functions on Q, with distributional 
first derivatives in L 2 (Q) provided with the norm 

IklU^tfi) = (livuiii^n) + hWhm) ■ (i-i) 

Here B is a ball at a positive distance from dil. If ft is unbounded, by L 1,2 ($7) we mean the completion of the 



space of functions with ||Vu||£2(Q) < oo, which have bounded supports, in the norm (1.1). The space of traces 
of functions in L l ' 2 {Vl) on dfl will be denoted by L 1 / 2 ' 2 (9ri). 

The maximum of diameters of F^\ j = 1, . . . , N is denoted by e. An array of points O^-*, j = 1, . . . , N, is 
chosen in such a way that O^' is an interior point of F^' for every j = 1, . . . , N. By 2d we denote the smallest 
distance between the points within the array {0"-'}v^_ 1 . It is assumed that there exists an open set lo C O 
situated at a positive distance from <9f2 and such that 

Uf =1 F {j) C L), dist( Uf =1 F {} \duj) > 2d, and diam u = 1, (1.2) 

With the last normalization of the size of ui, the parameters e and d can be considered as non-dimensional. The 
scaled open sets e~ 1 F^> are assumed to have Lipschitz boundaries, with Lipschitz characters independent of 
N. 

Our goal is to obtain an asymptotic approximation to a unique solution un € L 1,2 (Q,n) of the problem 

- A Ujv (x) = /(x) , xe^, (1.3) 

ujv(x)=^(x), xedft, (1.4) 

(x) = 0, xe5F«,j = l,...,7V, (1.5) 



9rj 



where <fi & i 1 / 2 ' 2 (951) and /(x) is a function in _L°°(f2) with compact support at a positive distance from the 
cloud u of small perforations. 

We need solutions to certain model problems in order to construct the approximation to ujy; these include 

1. v as the solution of the unperturbed problem in f2 (without voids), 

2. T>^ ""' as the vector function whose components are the dipole fields for the void F^ k \ 

3. H as the regular part of Green's function G in Vt. 

The approximation relies upon a certain algebraic system, incorporating the field Vf and integral character- 
istics associated with the small voids. We define 



0=1 * 

dxi 



«"")• £< 0,,, >' £< 0( ">- £' om '' £«"""'• fe<°"°>) T ■ 



and 6 = [S»j]i 7 -=i which is a 37V x 3iV matrix with 3x3 block entries 



®»i — "1 w=o tJ 



=0 w if * ^ J 



(V z <g>V w )(G(z,w)) 

O/3 otherwise 

where G is Green's function in $7, and -Z3 is the 3x3 identity matrix. We also use the block-diagonal matrix 

Q = diag{g«,...,gW}, (1.6) 

where Q} ' is the so-called 3x3 polarization matrix for the small void F^ k ' (see [7] and Appendix G of [8]). The 
shapes of the voids F»' , j = 1, . . . , N, are constrained in such a way that the maximal and minimal eigenvalues 
AmL, ^min of the matrices —Q^' satisfy the inequalities 

Aie 3 > max \$ ax , min X^l > A 2 e 3 , (1.7) 

l<j<N 1<]<N mm 



where A\ and Ai are positive and independent of e. 

One of the results, for the case when fl — M 3 , H = 0, and when (1.4 1 is replaced by the condition of decay 
of un at infinity, can be formulated as follows 

Theorem 1 Let 

e < cd , 

where c is a sufficiently small absolute constant. Then the solution ujy(x) admits the asymptotic representation 

N 

u N (x) = V (x) + J2 C(k) ■ ^ (fe) (x) + K N (x) , (1.8) 

fe=i 

where C (k) = {C[ k \ Cf \ cf } ) T and the column vector C = (cf \ cf \ cf\ . . . , c{ N) , C ( 2 N \ cf ] ) T satisfies 
the invertible linear algebraic system 

(I + 6Q)C = -0. (1.9) 

The remainder IZn satisfies the energy estimate 

\\VK N \\l 2(nN) < const {e 11 ^ 11 + e 5 d- 3 }\\Vv\\l Hn) . (1.10) 

We remark that since e and d are non-dimensional parameters, there is no dimensional mismatch in the right- 



hand side of (1.10 1. 



We now describe the plan of the article. In Section [21 we introduce the multiply-perforated geometry and 
consider the above model problems. The formal asymptotic algorithm for a cloud of small perforations in the 
infinite space and the analysis of the algebraic system ( |1.9[ ) are given in Sections [3] and |4j Section [5] presents the 
proof of Theorem [TJ The problem for a cloud of small perforations in a general domain is considered in Section 
[6j Finally, in Section [7] we give an illustrative example accompanied by the numerical simulation. 

2 Main notations and model boundary value problems 

Let n be a bounded domain in R 3 with a smooth boundary dil. We shall also consider the case when il = M 3 . 
The perforated domain H^v, is given by 



n N = n\uf =1 FU) , 

where F^' are small voids introduced in the previous section. Also in the previous section we introduced the 
notations e and d for two small parameters, characterizing the maximum of the diameters of F^> ,j = 1, . . . , N, 
and the minimal distance between the small voids, respectively. 

In sections where we are concerned with the energy estimates of the remainders produced by asymptotic 
approximations we frequently use the obvious estimate 

N < const gT 3 . (2.1) 

We consider the approximation of the function un which is a variational solution of the mixed problem 



(1.3H1.5J. 

Before constructing the approximation to ujv, we introduce model auxiliary functions which the asymptotic 
scheme relies upon. 

1. Solution v in the unperturbed domain fi. Let v G L l ' 2 {Vt) denote a unique variational solution of the 
problem 

-A«(x) =/(x), xeft, (2.2) 

u(x) = <^(x) , xean. (2.3) 



2. Regular part of Green's function in fl. By H we mean the regular part of Green's function G in O defined 
by the formula 

J ff(x,y) = (4^|x-y|)- 1 -G(x,y). (2.4) 

Then H is a variational solution of 

A x tf(x,y) = 0, xjell, 
tf(x,y) = (47r|x-y|)-\ xeffljefi. 

3. The dipole fields V { p , i = 1,2,3, associated with the void F^ . The vector functions T> u) = {X>- j) } 3 =1 , 
which are called the dipole fields, are variational solutions of the exterior Neumann problems 



A2> (j) (x) = 0, xeR 3 \F^, 

?-= — (x)=n«, xeSF®, 
on 

X>°' ) (x) = 0(e 3 |x-O(^r 2 



(2.5) 



where n") is the unit outward normal with respect to F^'. In the text below we also use the negative 
definite polarization matrix Q^' — {Qf k }? k=1 , as well as the following asymptotic result (see [7] and 
Appendix G in [8J, for every void F^>: 

Lemma 1 For |x— 0"'| > 2e, the dipole fields admit the asymptotic representation 

1 3 _ n^ 

^w4^ s ™^w + °( £4|x -° (,)| 1' l = w - (2 - 6) 

m— 1 ' ' 

The shapes of the voids F^> ,j = 1, . . . , TV, are constrained in such a way t hat the maximal and minimal 
eigenvalues Xmax, A^ n of the matrices — QS 3 ' satisfy the inequalities (1.7). 



3 The formal approximation of un for the infinite space containing 
many voids 

In this section we deduce formally the uniform asymptotic approximation of un'- 

N 

Ujv (x) ~ «(x) + Y, ° {k) ■ vW W . 
fc=l 

for the case Q = ]R 3 and derive an algebraic system for the coefficients C^ ' = {C\ } 3 =1 , k = 1, . . . ,N. 
The function un satisfies 

- Aujv(x) = /(x) , xefi„, (3.1) 



8un 



(x)=0, xeaF w ,j = l,...,JV, (3.2) 



07! 

uat(x) — > , as |x| — > oo . (3-3) 

We begin by constructing the asymptotic representation for u^ in this way 

N 

u;v(x) = v(x) + J2 Cik) ■ p(fc) ( x ) + K «W (3-4) 

where IZn is the remainder, and u(x) satisfies 

-Av(x) = /(x) , xel 3 , 



v(x) — > as |x| — > oo 



and V^ are the dipole fields defined as solutions of problems (2.5). The function TZn is harmonic in fi^r and 

K N {x) = 0(|x|- 1 ) as|x|^oo. (3.5) 



Placement of (3.4 1 into (3.2) together with (2.5) gives the boundary condition on dF^' 
dK N 



On 



-(x) = -n^-{Vv(0^) + C (j) +0(£)+ ^ V(C (fc) ■ 2? (fe) (x))}. 



Kk<N 



Now we use (2.6), for T>( k >, k ^ j, so that this boundary condition becomes 
d1l N 



On 



-(x) ~ -n^-{Vt>(0^) + C (j) + ^ r(x,0 (fe) )Q (fe) C (fc) }, xedF^\j = l,...,N , 



Kk<N 



where 



T(x,y) = (V z ®V w ) 



1 



^4ir\z — w 

Finally, Taylor's expansion of T(x, O^-*) about x = 0^\ j ^ fc, leads to 
dTZ N 



Z— X 

w=y 



(3.6) 



dn 



■(x)~-nW.{vi;(0 (j) ) + C (j) + ^ T(0«, 0^)Q w C (fc) } , xeaF«,j = l, 



,iV 



Kfc<Af 



To remove the leading order discrepancy in the above boundary condition, we require that the vector coefficients 
C"> satisfy the algebraic system 



Vu(O^) + C U) + J2 T(O w ,oW)Q (fc) C (fe) =0, ioij = l,...,N 



(3.7) 



Kk<N 



where the polarization matrices Qy' characterize the geometry of F^\ j = 1, . . . , N. Upon solving the above 
algebraic system, the for mal asymptotic approximation of u^ is complete. The next section addresses the 
solvability of the system (3.7), together with estimates for the vector coefficients C^. 



4 Algebraic system in the case Q = 1R 3 

The algebraic system for the coefficients C^ can be written in the form 

C+5QC= 0, 
where 

C = ((C«f, . . . , (cWff, © = ((v«(o«)f , . . . , (V,(0< s >)) T f, 

are vectors of the dimension 3iV, and 

1 



(4.1) 



& — [<%j Li = l> Sij — 



(V z <g>V w ) 



0/3 



47r|z — wl 



z=o« 

w=0< J > 



if i ^ j 



otherwise, 



Q = diag{Q (1) , . . . , Q (Af) } is negative definite. 



(4.2) 
(4.3) 



These are 2>N x 3N matrices whose entries are 3x3 blocks. The notation in (4.2) is interpreted as 

^ 3 



<->i 



f 1 d / Zr- ( r j) 

\4ird7 q \\z-0(ti\ 3 



K=OC« 



when i =/= j. 



q,r— 1 



We use the piecewise constant vector function 

f e a) C«), when x G B^ 4 , j = 1, ... ,7V, 
S(x) = I (4.4) 

l, 0, otherwise, 

where B { r j) = {x : |x - 0^| < r}. 

Theorem 2 Assume that X m ax < const d 3 , where A mQX is i/ie largest eigenvalue of the positive definite matrix 



— Q and f/ie constant is independent of d. Then the algebraic system (4.1 1 is solvable and the vector coefficients 
C^> satisfy the estimate 

^|( C 0)f Q«C«| < (1- const ^)- 2 ^|(V V (0^)) T g^Vt;(0«)|. (4.5) 



We consider the scalar product of (4.1 ) and the vector QC: 

(C, QC) + («SQC, QC) = -(©, QC). (4.6) 

Prior to the proof of Theorem we formulate and prove the following identity. 
Lemma 2 a) The scalar product («SQC, QC) admits the representation 

(SQCQC) =^[ I ^^(V.H(X))(V-H(Y))dYdX 

ir N 

3 ^|S^C«| 2 . (4.7) 



ird 3 



o,) TTie following estimate holds 



|(5QC,QC)| < const d~ 3 J^ IQ (j) C w | 2 , 

1<3<N 

where the constant in the right-hand side does not depend on d. 

Remark. Using the notation A/"(V • S) for the Newton's potential acting on V • S we can interpret the 



integral in (4.7) as 

(A/"(V-E),V-B 

V V h /L 2 (B3)' 

since obviously V • 3 <E W _1 ' 2 (R 3 ) and A/"(V • S) £ W 1,2 (M. 3 ). Here and in the sequel we use the notation (tp, ip) 
for the extension of the integral J m3 tp (X.)ip(K)dX onto the Cartesian product W 1:2 (M 3 ) x W~ 1,2 (R 3 ). 
Proof of Lemma [2l a) By ( |4.2[ ), ( |4.3| , the following representation holds 

(SQC,QC) = lf:(Q^) T £ (V^V w )(^)|^ o( ,(Q( fe ) C W). (4.8) 

j = l l<k<N,kjtj Vl l/ w=0 (fc) 

Using the mean value theorem for harmonic functions we note that when j 7^ k 



1 



(V Z ®V W ) 1 r 0) = / (V z <g>V^ 



1 



dw. 



Substituting this identity into (4.8 1 and using definition (4.4) we see that the inner sum on the right-hand side 



of (4.8) can be presented in the form 

Trd 3 T^S+J V \ B w idY q \\ 



- fejL fl «> {^(|Y-5)|»)}lr=l S(Y)dY ' 



(i 



and further integration by parts gives 



(5QC,QC) 



3=1 

{[ ( Y — O® ~i 3 



(4.9) 



dY 



■>(j) 



(Jh 



|Y-OW)|=(d/4)— 7 



(Y r -Q¥>)(Y q -O q J >) 
|Y-O0')|4 



dSy Q U) C U) 



r,q=l 



where the integral over K J \ B(d/4)_ T in (4.9) is understood in the sense of distributions. The surface integral 
in (4.9) can be evaluated explicitly, i.e. 



|Y-OW)|=(d/4)-r 



(Y r - Q^)(Y q - O q j) ) ' 
|Y-0«l 4 



r,q— 1 






(4.10) 



Once again, applying the mean value theorem for harmonic functions in the outer sum of (4.9) and using (4.10) 



together with the definition (4.4) we arrive at 

N 



(5QC,QC) 



1G 



-^ KmVf / Vs r (X)-^f , * , W-B(Y)dYdX, 



(4.11) 



j = l"' S (d/4) + r ,;K3 \- B (d/4)-T r=l 



where H r are the components of the vector function S defined in (4.4) 



The last integral is understood in the sense of distributions. Referring to the definition (4.4), integrating by 

(4.12) 



parts, and taking the limit as r — > 0+ we deduce that the integral term in (4.11) can be written as 

576 



7T 3 d 6 



1 



IY-XI 



V-S(X) V-S(Y) )dYdX 



Using (4.11) and (4.12) we arrive at (4.7). 



b) Let us introduce a piece-wise constant function 



c(x) = J cU) ' when x G B dh ' i = !. 
, otherwise . 



,JV 



According to the system (3.7), V x C(x) = O, and one can use the representation 

C(x) = VW(x) 



(4.13) 



where W is a scalar function with compact support, and (4.13) is understood in the sense of distributions. We 

-1 
4- 



give a proof for the case when all voids are spherical, of diameter e, and hence Q 3 ' = — j£ 3 l3, where ^3 is the 



identity matrix. Then according to (4.11) we have 



|(5QC,QC)| < ^ i EjQ°" ) ^ ) | 2 + ^|/ R3 / K3 (Vx^ ( X).V x (^))A Y M,(Y ) ^X 

J2 I U) |VT^(Y)| 2 dY 



l<j<N 



^EiW + ¥ 



< 



7TG? 3 

l<j<N 

const 



l<j<N " B d/i 



d 3 



£ \Q U) c 



0)1 2 



l<j<N 



□ 

Proof of Theorem [2l Consider the equation (4.6). The absolute value of its right-hand side does not 



exceed 

(C-QC^e.-Qe) 1 ' 2 . 

Using Lemma 1 and part b) of Lemma 2 we derive 

(C, -QC) - const^ 3 (-QC, -QC) < (C, -QC)^ 2 (&, -Q0) 1/2 , 
leading to 

( 1 -^t§^Qp )( c,_ QC) .^ (9 ,_ Q e,v», 

which implies 

(l- const ^^) 2 (C,-QC) < (0, Q0). (4.14) 

The proof is complete. □ 

Assuming that the eigenvalues of the matrices —Q} J ' are strictly positive and satisfy the inequality (1.7), 
we also find that Theorem [2] yields 



Corollary 1 Assume that the inequalities (1.7) hold for X max and X m in- Then the vector coefficients C^' i 



in 



the system (4.1) satisfy the estimate 



]T |C«| 2 < const cr 3 ||V V ||i 2(w) , (4.15) 

l<j<N 



where the constant depends only on the coefficients A\ and A 2 in (1.7) 



Proof. According to the inequality (4.5) of Theorem [2] we deduce 



X min J2 |C°Y < (1 - ^A™*)- 2 ^ J2 |V«(0«)| 2 . (4.16) 

l<j<N l<j<JV 

We note that v is harmonic in a neighbourhood of to. Applying the mean value theorem for harmonic functions 
together with the Cauchy inequality we write 

\ Vv{0 U)^< || Vu || 2 

7rcr L 2\ B d/4> 



Hence, it follows from (4.16) that 



y: ic«i 2 <d- 3 (i-^A mo ,r 2 ^^ x: nv<« W) 

l<j<N mm l<j<N ' 

SA-*(<-\ COnst \ \-2 48 \noaA ,,„ , |2 ,„ -, 7 n 

<d 1(1 ^-A max ) — ||Vv|| i2M , (4.17) 

\ w 7T strain. I 



which is the required estimate (4.15). D 



5 Energy error estimate in the case £1 = IR 3 

In this section we prove the result concerning the asymptotic approximation of ujy for the perforated domain 
51 at = R 3 \U^_ 1 .FC'). The changes in the argument, necessary for the treatment of a general domain, will be 
described in Section |6] 



Proof of Theoremni a) Neumann problem for the remainder. The remainder term 72 at in ( 1.8 ) is a harmonic 
function in Qn, which vanishes at infinity and satisfies the boundary conditions 



<972 



8 



^(x) = -(Vi;(x)H-C W ") -n«(x)- V C (fe) - ilz>(*)(x), when x e dF&, j = 1, . . . ,N. (5.1) 
an V / f— J on 



k+i 
Kk<N 



Since supp / is separated from F^\j = 1, . . . , N, and since T>^\ j = 1, . . . , N, satisfy (2.5) we have 

dn N 



8F 



U) dn 



(x)dS x = 0, j = l,.,.,N, 



(5.2) 



b) Auxiliary functions. Throughout the proof we use the notation B p = {x : |x — O^'l < p}. We introduce 
auxiliary functions which will help us to obtain (1.10). Let 



tffc(x) = I i(x)-«(OW)-(x-OC ; )).V«(OW)+ Y C (j) -X» (j) (x) 



l<j<N 



Y (x -O^) -T(0 (fc) ,O w )Q U) C U) , 



(5.3) 



l<j<N 
3^k 



for all x <E 0^ and fc = 1, . . . , JV. Every function $j. satisfies 

-A^ fe (x) = /(x), xen Nl (5.4) 

and since u n supp / = 0, we see that \tfc, fe = 1, . . . , N, are harmonic in to. Since the coefficients C^ satisfy 



system (4.1), we obtain 



£« + ^M-o, ««*■> 



and according to (5.2) the functions &}. have zero flux through the boundaries of small voids F^ k \ i.e. 

f d^k . 

J dF (k) dn 

Next, we introduce smooth cutoff functions 



-(x)dx = ,k = l,...,N 



(5.5) 



(5.6) 



,(*>) 



x-> x ((x-0«)/ £ ), k = l,...,N, 



equal to 1 on B 2e and vanishing outside B^J . Then by (|5.5|) we have 



d_ 

On 



(n N (x)+ J2 xi fe) (x)* fc (x))=0 on dF®, j = l,. 



,N. 



(5.7) 



Kk<N 



c) Estimate of the energy integral oflZjy in terms of^k- Integrating by parts in $7^ and using the definition 
of Xs > we write the identity 

f VK n -v(k n + Y X { e k) *k)d X = - f TZ N A(n N + Y xi k) *k)d X , (5.8) 

which is equivalent to 

/ \VK N \ 2 d x + Y, I fcl _ V^-V(x«* fc )dx = - E / rfel _ m ^A( X W* fc )dx, (5.9) 



Kk<N" a 3 



Kk<N" a 3 



since 7?. at is harmonic in fijv. 

We preserve the notation 72 n for an extension of IZn onto the union of voids F^ with preservation of the 

class W 1 ' 2 . Such an extension can be constructed by using only values of 72. jv ° n the sets B^J \F in such a 
way that 



IV72 



N\\ L2[B W 



< const V72 



N "L2(B£\F (k) )- 



(5.10) 



The above fact follows by dilation x — > x/e from the well-known extension theorem for domains with Lipschitz 
boundaries (see Section 3 of Chapter 6 in [10]). We shall use the notation 1Z for the mean value of 1Z n on 



B 



(fe) 



3e 



The integral on the right-hand side of (5.9 1 can be written as 



Y f Tl N A( X i k ^ k )dx=- y [ (K N -K {k) )A( 

ik<N Jts -iE \ b Kk<N Jts 3c \* 



{k) ^( X ^ k )d^ 



In the derivation of (|5.11|) we have used that 

A 

B W\pW 



(x^ k ) 



dx — 



i)F 



cm 

(k) dn 



dS* = 



(5.11) 



(5.12) 



,(*0 



Owing to (5.8) and (5.11), we can write 



according to (5.6) and the definition of \ 

te 

l|V^Ar||i 2(njv) <S 1+ S 2 , 



where 



Kk<N "^3 



and 



s 2 = y I / (K N -n {k) )A( x i k \* k -* k ))dx 



Kk<N " D 3 



(5.13) 
(5.14) 

(5.15) 



>(*) 



where Vfrfc is the mean value of Vffc over the ball B\ . Here, we have taken into account that by harmonicity of 



,(.*>) 



IZn, (|5.2|) and definition of xl 

f A(ll N -ll {k) )xi k Ux= J A(K N -K (k) )xi k) dx = 0. 

J b v>\fM V J J B w V J 



By the Cauchy inequality, the first sum in (5.13) allows for the estimate 



Ei<( y IIVT^II^^, 



l<fe<AT 

Furthermore, using the inequality 



1/2 



( E H 



xi fe) * fc 



Kk<N 



L 2 {B < £ > \F W \ 



1/2 



XI HV^H^fB^N^e^ ^ ll V72 -iv|||2 (njv) , 



Kk<N 



together with (5.16), we deduce 



^i<\\VTZ N \\LHn N )( y ||V(xW**j||,„, .,,, r .,., 



Kfc<JV 



2 N ! -' 

i2 ( B 3e ; \-P""')> 



(5.16) 



(5.17) 



(5.18) 



Similarly to (5.16), the second sum in (5.13) can be estimated as 



S 2 < 



E 



Kk<N " a 3e 



i (n N -n (k) fd x ) 1/2 ( f (A (x «(^-^.))) 2 d x) 1/2 . (5.19) 



By the Poincare inequality for the ball B 



(k) 

3e 



l|ft"-tt ( % (flS )) 



B W\ F v 



<conste 2 ||V^Ar|| 2 r2 , D(fc)x 



'£»(B^) 



(5.20) 
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we obtain 



E 2 < const e ( ^2 ||V7e.-- :; 

Kk<N 



N »L*(B^) 



1/2 



E 



l<k<N" B 3e \* 



(A(xi fc) (vl/ fe -vl/ fc ))) 2 dx 



which does not exceed 



const £ ||V^|U 2(ajv) ( X! /,„ _ (fc) (A(xi fc) (* fc -*fc))) 2 ^ 



1/2 



(5.21) 



because of (5.101. Combining ( 5.13[ )-( |5~2"l ) and dividing both sides of (5.13) by \\VR-n\\l 2 (q n ) we arrive at 

1/2 



\VTl N \\ LHnN) <( J2 V(x^(*fc-*fc)) 

Kk<N 



Z/"(u£>) 



+const£(^ /" (fc) {(* fe -* fc )Axi' £) +2Vxi fe) -V* fc } 2 dx) 1/2 , (5.22) 



KKJV" D 3e 



which leads to 



ivk 



JVllz^fijv) 



< const 



E (ll V ^l| 2 L 2 (^))+ e " 2 H^-^llL 2 (B«) 



Kfe<iV 



(AO 



Applying the Poincare inequality (see (5.20)) for "J^. in the ball B\ e and using (5.23), we deduce 



l|V^Ar||| 2(0jv) < const J2 H V * 



Kk<N 



k|l L 2( B W 



(5.23) 



(5.24) 



rfj Final energy estimate. Here we prove the inequality (1.10). Using definition (5.3) of ^k, k = 1, . . . ,N, 
we can replace the preceding inequality by 



|V^7v|||2 (njv) < const {/C + £}, 



(5.25) 



where 



K= J2 ||Vu(-)-VKO (fc) )|| 



Kk<N 



L*(B<*>) ' 



£ = 



E ii E [V(C« •©«(•)) -T(0 (fe ),0«)Q w C« 



l</c<AT j^fc 
l<j<N 



(5.26) 



l»(b£Y 



The estimate for K. is straightforward and it follows by Taylor's expansions of v in the vicinity of 0^ k \ 

d 2 v 2 



/C < const e 5 d 3 _max 

xGw,l<i,j<3 



dXidXj 



(5.27) 



Since w is harmonic in a neighbourhood of w, we obtain by the local regularity property of harmonic functions 
that 

(5.28) 



K < const e 5 d~ 3 \\Vv 



m 



To estimate £, we use Lemma [lion the asymptotics of the dipole fields together with the definition (3.6) of 
the matrix function T, which lead to 



|V(C W -£>W(x)) -T(0 (fc) ,0 (j '))Q (i) C w | < const £ 4 |C w ||x - O 



0)1-4 



>(*0 



for x e B^. Now, it follows from (5.26) and (5.29) that 



£ < const £ 8 V / f V , ' -tL ) dx, 



(5.29) 



(5.30) 
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and by the Cauchy inequality the right-hand side does not exceed 

N N N 



N N r A N N 

const £ 8 V|CW| 2 y V / **— -< const e 11 V |C^| 2 V V , . n * ...,„ 

p=l k=ll<j<N,j^k J±l ^ ' ' p=l fe=l l<i<JV,j'^fe I I 



-11 



s »- £ £ ^Jl,„ di ct * — V £ icH 



p=i 



11 w 



p=l 



Since the eigenvalues of the matrix — Q satisfy the constraint (1.7 1, we can apply Corollary [T] and use the 
estimate (4.15) for the right-hand side of (5.31) to obtain 



£< const e ll d~ 11 II Vdl? 



L 2 (u 



(5.32) 



Combining (5.25), (5.28) and (5.32), we arrive at (1.10) and complete the proof. □ 



6 Approximation of un for a perforated domain 



Now we seek an approximation of the solution ujy to the problem (jl .31) — (1.5) assuming that Q, is an arbitrary 



domain in R . We first describe the formal asymptotic algorithm and derive a system of algebraic equations, 



similar to (4.1), which is used for evaluation of the coefficients in the asymptotic representation of um- 



6.1 Formal asymptotic algorithm for the perforated domain Q^ 

The solution un G L 1,2 (CIn) of (1.3)— (1.5) is sought in the form 



UN 



N 

(x) = «(x)+EC W ' {£> (fc) M - Q (fe) V y #(x,y)| y=0(fc) } + R N (x) 



(6.1) 



fc=i 



where in this instance v solves problem (2.2), (2.3) in Sectional and Rn is a harmonic function in f2jv- Here 
C' ' , k = 1, . . . , N are the vector coefficients to be determined. 

Owing to the definitions of T>^ ', k — 1,...,N, and H as solutions of Problems 2 and 3 in Section [2j and 
taking into account Lemmalllon the asymptotics of T >^ ' w e de duce that |i?jv(x)| is small for x G dfl. 

On the boundaries dF^\ the substitution of (6.1) into (1.5) yields 



dR 



N 



On 



(x) = -n w ■ hv(0^) + C U) + 0(e) + 0{e 3 \C U) \) 

+ J2 v{cW-(l?W(x)-Q( fe )v y ff(x, y )| y=0(fc) )}}, x£Sf«i = l,..., 



N . 



Kk<N 



Then, using the asymptotic representation (2.6) in Lemma [T] we deduce 
dR N 



dn 



■(x) 



,(;;) 



{Vu(0 (i) ) + C W + J^ 3:(x,0( fc) )S (fc) C (fc) } , xG9F«,j = l,...,iV, (6.2) 



Kk<N 



(6.3) 



where T(x, y) is defined by 

T(x,y) = (V x ®V y )G(x,y), 

with G(x, y) being Green's function for t he d omain fi, as defined in Section [2] To compensate for the leading 
discrepancy in the boundary conditions ( |6.2| , we choose the coefficients C'^*, m — 1, . . . , N, subject to the 
algebraic system 

Vv(0^) + C (j) + J^ ^(0 U) ,0 {k) )Q {k) C (k) =0, j = l,...,N, (6.4) 



Kk<N 



where Q^ ', k = 1, . . . , N, are polarization matrices of small voids F^ , as in Lemma [lj 
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Provided system (6.4) has been solved for the vector coefficients C^- k \ formula (6.1) leads to the formal 

(6.5) 



A 



asymptotic approximation of un'- 

u N (x) ~ „(x) + J2C W ■ {^ (fe) (x) - Q (fc) V y i/(x,y)| y=0(fe) } 

6.2 Algebraic system 



fc=i 



The system (6.4) can be written in the matrix form 

C + 6QC 
where 



-e, 



© = [6y]£ i=1 , © 



(V z ®V w )G(z,w) 



Oh 



_ w if * ^ 3 



otherwise 



(6.6) 



(6.7) 



with G(z,w) standing for Green's function in the limit domain fi, and the block-diagonal matrix Q being the 
same as in (1.6). The system (6.6) is similar to that in SectionBJ with the only change of the matrix S for ©. 
The elements of © are given via the second-order derivatives of Green's function in 51, as defined in (6.3). The 
next assertion is similar to Corollary [T] 

Lemma 3 Assume that inequalities ( |1.7[ ) hold for X max and X m in- Also let v be a unique solution of problem 
(2.2), (2.3) in the domain O. Then the vector coefficients C^' in the system (6.4) satisfy the estimate 



J2 \C U) \ 2 < const cT 3 ||Vv|| 

l<j<JV 



2 

L 2 (n)' 



(6.8) 



where the constant depends on the shape of the voids F^',j — 1, . . . , N. 



Proof. The proof of the theorem is very similar to the one given in Section [4] We consider the scalar product 
of (6.6) and the vector QC: 

(C,QC) + (©QC,QC) = -(0,QC), (6.9) 



and similarly to (4.7) derive 



(6QC,QC) =48 2 tt- 2 (T 6 / / G(X, Y)(V ■ S(X))(V • S(Y))rfYdX 

Jn Jn 

-Wir^d- 3 Yl IQ (i) C (i) | 2 

l<j<JV 

J2 (q^cw) T (v z ®v w )(^(z,w))| z=o0) (q«c«), 



l<j<N 



(6.10) 



-oO) 



w=0 



where the integral in the right-hand side is positive, and it is understood in the sense of distributions, in the 
same way as in the proof of Lemma ^1 while the magnitude of the last sum in (6.10) is small compared to the 
magnitude of the second sum. 

Now, the right-hand side in (|6.9|) does not exceed 



(C,-QC) 1 / 2 (&,-Q&) 1 / 2 . 
Following the same pattern as in the proof of Theorem [2] we deduce 

(C, -QC) - const d- 3 (-QC, -QC) < (C, -QC)^ 2 (&, -Q0) 1 / 2 
where the constant is independent of d. Furthermore, this leads to 

i 1 C ° nSt d ' 3 ( '(C C -QC) C) ) <C ' - QC>1/2 " ( ®' - Q0>1/2 ' 
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which implies 



(l - const <T 3 A roc « ) (C. - -QC) _. (0. - QO; 



(6.11) 



where \ max is the largest eigenvalue of the positive definite matrix — Q. Then using the same estimates (4.16) 
and (4.17) as in the proof of Corollary 111 we arrive at (6.8). □ 



6.3 Energy estimate for the remainder 
Theorem 3 Let the 'parameters e and d satisfy the inequality 

e < cd . 



where c is a sufficiently small absolute constant. Then the solution ujy(x) of ( 1.3 )— ( 1.5 ) is represented by the 
asymptotic formula 



N 



U N 



(x) = «(x) + £ C« • {2>«(x) - Q™ V y H(x, y) } + R N (x) , 



fc=i 



y=0( fc )' 



(6.12) 



where C^ ' = (Cj , C\ \C\ ) solve the linear algebraic system (6.4). The remainder Rn in (6.12) satisfies 
the energy estimate 

(6.13) 



\VR N \\ 2 L2{nN) < const { £ 11 d- 11 + £ 5 d- 3 }||V^| 



2 

L 2 (n) • 



Proof. Essentially, the proof follows the same steps as in Theorem [II Thus, we give an outline indicating 
the obvious modifications, which are brought by the boundary dil. 

a) Auxiliary functions. Let us preserve the notations Xs for cutoff functions used in the proof of Theorem 
[l] We also need a new cutoff function xo to isolate dQ from the cloud of holes. Namely, let (1 — Xo) S Co°(^) 
and xo = on a neighbourhood of U. A neighbourhood of dft containing supp xo will be denoted by V. Instead 



of the functions S£>k defined in (5.3), we introduce 

*i n) (x) = w(x)-w(0 (fc) )-(x-0 ( ' £) )-Vw(0 (fc) )+ J2 C U) -T> U \x) 



l<j<N 

N 



J2 (x- O w ) • I(0 ( * ) ,0 W) )fi w C (,-) - J]C (J) • Q (j) V y ff(x,y) 

3=1 



l<j<N 



=o«> 



(6.14) 



and the algebraic system (6.4) we have 



where the matrix X is defined in (6.3) via second-order derivatives of Green's function in fi. Owing to (6.12) 



0_ 

On 



(*£ n) (x) + fljv(x)) = 0, x G 9F (fc) . 



We also use the function 



N 



i=i ' 



(6.15) 



(6.16) 



which is harmonic in fijy. It follows from (6.12) that 



(x - Ok')) 



JMx) + *o(*) = - E Cb) • QiJ) LS~- OW - Vyg(X ' Y t=o J = ° ' xe9 " 



l<j<N 



(6.17) 
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b) The energy estimate for Rn ■ We start with the identity 

f v(r n + X o%)-v(Rn+ Y, xi fc) *i n) ) dx 

J&N l<k<N 

= -[ (r n + Xo^o)a[r n + J2 X$! l) *$ 1) )d*, 



(6.18) 



Kk<N 



ik) 



which follows from (6.15), (6.171 by Green's formula. According to the d efinitions of xo and Xe ■> we have 
supp xo H supp Xe = for all fc = 1, . . . , N. Hence the integrals in (6.18) involving th e pro ducts of xo and 

(k) J I 

Xe or their derivatives are equal to zero. Thus, using that ARn = on J7jv, we reduce (6.18) to the equality 



/ |Vi?Ar| 2 dx+ V f Vi? N -v(xi fc) *i a) )dx+ f VR N -vfxo^oW (6.19) 



E 



Kk<N " ' B 3e\F 



R N A( x i k ^)dx, 



Kk<N" ilN 



which differs in the left-hand side from (5.9) only by the integral over f2jv H V. 
Similarly to the part (b) of the proof of Theorem IT] we deduce 



|V^Ar||| 2(njv) <const{||V* |l! 2 



(Onv) + ll*o||.L 2 (nnv) 



E ii v ^iii= ( B«)} 



Kk<N 



Similar to the steps of part (d) of the proof in Theorem n] the last sum is majorized by 

const (e n rf- n +e rj d- 3 )\\Vv\\ 2 L 2 {n) . 



It remains to estimate two terms in (6.20) containing ty . Using ( 5.29[ ), together with (6.8) we deduce 

mt ..a s V- f \C {J) \ 2 < 



i<i<JV' 



\C^\ 2 dx 

)0')|6 



< const e 8 J2 \ CU) \ 2 ^ const^||Vw| 
i<i<N 



2 

L 2 (n)> 



and 



l|V*o||i»(OnV) < const e 8 £ / -J^J- 

^.<- H ^nv |x-u 



1<J<JV 



|C^)| 2 dx 

0)|8 



< const e 8 J2 IC '^ 2 < const— ||Vv| 



2 
L 2 (f2)- 



Combining ( 6.20 H 6.23 1 we complete the proof. □ 



(6.20) 



(6.21) 



(6.22) 



(6.23) 



7 Illustrative example 

Now, the asymptotic approximation derived in the previous section is applied to the case of a relatively simple 



geometry, where all the terms in the formula (6.12) can be written explicitly. 



7.1 The case of a domain with a cloud of spherical voids 

Let H,n be a ball of a finite radius R, with the centre at the origin, containing N spherical voids F^' of radii 
Pj with the centres at 0^',j = 1, . . . ,N, as shown in Fig. Ill The radii of the voids are assumed to be smaller 
than the distance between nearest neighbours. We put = and 



/(x) 



6 when |x| < p, 
when p < |x| < R. 



(7.1) 
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Figure 1: Example configuration of a sphere containing a cloud of spherical voids in a the cube ui. 



Here, it is assumed that p + b < \0^\ < R— b, 1 < j < N, where p and b are positive constants independent 
of e and d. 

The function un is the solution of the mixed boundary value problem for the Poisson equation: 



Amat(x) + /(x) = 0, when x £ fl^, 
un(x) — 0, when |x| = R, 

^(x) = 0, when |x - 0&\ = Pj , j = l,.. 



,N. 



(7.2) 
(7.3) 

(7.4) 



In this case, un is approximated by ( |6.12[ ), where the solution of the Dirichlet problem in $7 is given by 

„^^ - f p 2 (3-2pi?- 1 ) - |x| 2 when|x|<p, 

W W _ 1 2p 3 (|x|- 1 - i?- 1 ) when p < |x| < R. [ > 



In turn, the dipole fields T>^' and the dipole matrices Qy' have the form 

x - OW 



X> (j) (x) = 



J |X-0(J)| 3: 



Q 



(j) _ 



-2np 3 J 3 , 



(7.6) 



where ^3 is the 3x3 identity matrix. 



The regular part H(x,y) of Green's function in the domain fl (see (2.4)) is 

R , R 2 



#(*,y) 



47r|y||x-y| 



y |yP y ' 



(7.7) 



The coefficients C^\ j = 1, 
G(x, y) is given by 



, N, in (6.12 ) are defined from the algebraic system (6.4), where Green's function 
f 



G(x,y) 



R 



47r|x-y| 47r|y||x-y| 



(7.8) 



7.2 Finite elements simulation versus the asymptotic approximation 

The explicit representations of the fields v,T>^' ,H,G, given above, are used in the asymptotic formula (6.12). 



Here, we present a comparison between the results of an independent Finite Element computation, produced in 



COMSOL, and the mesoscale asymptotic approximation (6.12) 



For the computational example, we set R — 120, and consider a cloud of N — 18 spherical voids arranged 
into a cloud of a parallelipiped shape. The position of the centre and radius of each void is included in Table 
fl] The support of the function / (sec (7.1 1), is chosen to be inside the sphere with radius p = 30 and centre at 
the origin, as stated in (7.1). 
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Void 


Centre 


Pi/R 


Void 


Centre 


Pi/R 


F (D 


(-50, 0, 0) 


0.0417 


p{W) 


(-72, 0, 0) 


0.0417 


i?(2) 


(-50, 0, 22) 


0.0333 


i?(n) 


(-72, 0, 22) 


0.0458 


p(3) 


(-50, 22, 0) 


0.0292 


i^(12) 


(-72, 22, 0) 


0.0292 


pW 


(-50, 0, -22) 


0.0375 


i?(13J 


(-72, 0, -22) 


0.0375 


i?(5) 


(-50, -22, 0) 


0.0458 


i?(14) 


(-72, -22, 0) 


0.0417 


i?(6) 


(-50, 22, 22) 


0.0292 


i?(15) 


(-72, 22, 22) 


0.0333 


i?(7) 


(-50, 22, -22) 


0.025 


^(16) 


(-72, 22, -22) 


0.05 


i?(8) 


(-50, -22, 22) 


0.0375 


i?(l7) 


(-72, -22, 22) 


0.0333 


i?(9) 


(-50, -22, -22) 


0.0375 


i?(18) 


(-72, -22, -22) 


0.0375 



Table 1: Data for the voids F<J\ j = 1, . . . , 18. 



Figure shows the asymptotic solution um of the mixed boundary value problem (part (b) of the figure) 
and its numerical counterpart obtained in COMSOL 3.5 (part (a) of the figure). This computation has been 
produced for a spherical body containing 18 small voids defined in Table [T] The relative error for the chosen 
configuration does not exceed 2%, which confirms a very good agreement between the asymptotic and numerical 
results, which are visually indistinguishable in Fig. [2^, and Fig. [2b. 




100 \/ 



r 



™0 




(a) 



(b) 



Figure 2: Perforated domain containing 18 holes: (a) Numerical solutions produced in COMSOL; (b) Asymptotic 
approximation. 

The computation was performed on Apple Mac, with 4Gb of RAM, and the number N = 18 was chosen 
because any further increase in the number of voids resulted in a large three-dimensional computation, which 
exceeded the amount of available memory. Although, increase in RAM can allow for a larger computation, it is 
evident that three-dimensional finite element computations for a mesoscale geometry have serious limitations. 
On the other hand, the analytical asymptotic formula can still be used on the same computer for sigtnificantly 
larger number of voids. 

In the next subsection, we show such an example where the number of voids within the mesoscale cloud runs 
upto N = 1000, which would simply be unachievable in a finite element computation in COMSOL 3.5 with the 
same amount of RAM available. 



7.3 Non-uniform cloud containing a large number of spherical voids 



Here we consider the same mixed boundary value problem as in Section |7.1| but the cloud of voids is chosen 
in such a way that the number N may be large and voids of different radii are distributed in a non-uniform 
arrangement. For different values of N, the overall volume of voids is preserved - examples of the clouds used 
here are shown in Fig. [T] 



The results are based on the numerical implementation of formula (6.12| in MATLAB. 



The cloud to is assumed to be the cube with side length -4= and the centre at (3,0,0). Positioning of voids 
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is described as follows. Assume we have N — m s voids, where m — 2,3, ... . Then u is divided into N smaller 
cubes of side length h = —i— , and the centres of voids are placed at 

° V3771 



Q (j>,q,r) = f 3 



1 



2p- 1 



h, - 



1 



2q-l 



2r-l 



2y/3 2 2^3 

for p,q,r — 1, . . . , m, and we assign their radii p p , q , r by 

h 
5, 



2^3 



Pp , g , r 



= < 



if p > q 
if p < q 
iip = q 



where a < 1, and it is chosen in such a way that the overall volume of all voids within the cloud remains 
constant for different N. An elementary calculation suggests that there will be m 2 voids with radius \ and 

equal number m ^~ ' of voids with radius § or ~ . 

Assuming that the volume fraction of all voids within the cube is equal to j3, we have 



4tt/i 3 /ra 2 (m - 1)(8 + 125a 3 ) 
3 V 



2000 



m 

G4 



= /3 



3^3' 



and hence 



In particular, if N 



a 



16m 
m 

oo, the limit value a„a becomes 






«; 



L 7T 



125 + 32(m - 1) 
8000m . 



8 -l 1/3 
125/ 



(7.9) 



(7.10) 



In the numerical computation of this section, j3 = 7r/25. 

Takin g R = 7 and p = 2, we compute the leading order approximation of uj^ — v, as defined in the asymptotic 
formula (6.12 1, along the line 7 at the intersection of the planes xi = — 1/(2\/3) and x% = — 1/(2\/3), for 
N = 8, 125, 1000. Fig. [3] below shows the configuration of the cloud of voids for a) N — 8 and b) N — 125. 
For a large number of voids (N = 1000), Fig. Bk) shows the cloud and FigEb) includes the graph of a versus 
N. The plot of un — v given by (6.12) for 2 < x\ < 4 is shown in Fig. [5J The asymptotic correction has been 
computed along the straight line 7 = {x\ € R, x 2 = — l/(2-v/3),x 3 = — l/(2y3)}. Dipole type fluctuations are 
clearly visible on the diagram. Beyond N — 1000 the graphs are visually indistinguishable and hence the values 
N = 8, 125, 1000, as in Figures [3] and [4] have been chosen in the computations. The algorithm is fast and does 
not impose periodicity constraints on the array of small voids. 

Acknowledgments. We would like to acknowledge the financial support of the U.K. Engineering and Physical 
Sciences Research Council through the research grant EP/F005563/1. 
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Cloud of voids. N=8 



Cloud of voids, N-125 





a) b) 

Figure 3: The cloud of voids for the cases when a) N — 8 and b) N = 125. 



Cloud of voids. N=1000 





000 2000 3000 4000 5000 6000 7000 8000 



a) b) 

Figure 4: a) The cloud of voids for the cases when N = 1000, b) The graph of a versus N given by formula 



(7.9) when (3 = it/25, for large N we see that a tends to 0.7465 which is predicted value present in (7.10) 
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2.2 2.4 2.6 2.8 



Figure 5: The graph of un — v given by (6.12 1, for 2 < x\ < 4 plotted along the straight line 7 adjacent to the 
cloud of small voids. 
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